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A general approach for calculating spectral and optical properties of pigment-protein complexes of known 
atomic structure is presented. The method, that combines molecular dynamics simulations, quantum chemistry 
calculations and statistical mechanical modeling, is demonstrated by calculating the absorption and circular 
dichroism spectra of the B800-B850 bacteriochlorophylls of the LH2 antenna complex from Rs. molischianum 
at room temperature. The calculated spectra are found to be in good agreement with the available experimental 
results. The calculations reveal that the broadening of the B800 band is mainly caused by the interactions with 
the polar protein environment, while the broadening of the B850 band is due to the excitonic interactions. Since 
it contains no fitting parameters, in principle, the proposed method can be used to predict optical spectra of 
arbitrary pigment-protein complexes of known structure. 



I. INTRODUCTION 

Pigment-protein complexes (PPCs) play an important role 
in photosynthetically active biological systems and have been 
the subject of numerous experimental and theoretical studies ^ . 
In a PPC the photoactive pigment molecules are held in well 
defined spatial configuration and orientation by a scaffold of 
proteins. The availability of high resolution crystal struc- 
ture for a continuously growing number of PPCs provides 
a unique opportunity in better understanding their properties 
and function at atomic level. The spectral and optical prop- 
erties of PPCs are determined by (i) the chemical nature of 
the pigment, (ii) the electronic interactions between the pig- 
ment molecules, and (iii) the interactions between pigment 
molecules and their environment (e.g., protein, lipid, solvent). 
Since in biological systems PPCs exist and function at physio- 
logical temperature their electronic and optical properties are 
strongly affected by thermal fluctuations which represent the 
main source of dynamic disorder in these systems. 

Unfortunately, even in the simplest theoretical models of 
PPCs the simultaneous treatment of the electronic coupling 
between the pigments and the effect of thermal disorder can 
be done only approximately^'^ The purpose of this paper 
is to formulate and implement an efficient method for calcu- 
lating spectral and optical properties, e.g., linear absorption 
(OD) and circular dichroism (CD) spectra of PPCs at finite 
temperature by using only atomic structure information. To 
demonstrate its usefulness, we apply the proposed method to 
calculate the OD and CD spectra at room temperature of the 
aggregate of bacteriochlorophyll-a (BChl-a) molecules in the 
LH2 antenna complex from the purple bacterium Rs. molischi- 
anum. Following their crystal structure determination, LH2 
complexes from molischianum^ and Rsp. acidophilcAhave 
been extensively studied both experimentalljs^iSiiSiiiiiSiiiiiii^ 
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and dieoretically^i^i^°i'^"i^^i'^i^°i^^^^i^^ InRs. molischianum 
the LH2 is an octamer of aj3-heterodimers arranged in a ring- 
like structure^'*'^^. Each protomer consists of an a- and a )3- 
apoprotein which binds non-covalently one BChl-a molecule 
that absorbs at 800 nm (refeiTed to as B800), two BChl-a 
molecules that absorb at 850 nm (referred to as B850) and 
at least one carotenoid that absorbs around 500 nm. The 
total of 16 B850 and 8 B800 BChls form two circular ag- 
gregates, both oriented parallel to the surface of the mem- 
brane. The excitonic coupling between the B800s is negligi- 
ble because of their large spatial separation (^ 22 A). There- 
fore, the optically active Qy excited electronic states of the 
B800s are almost degenerate. On the other hand, the tightly 
packed B850s (with and average Mg— Mg distance of ^ 9.2 A 
within the aj3— heterodimer and ~ 8.9 A between the neigh- 
boring protomers) are strongly coupled and the corresponding 
Qy excited states form an excitonic band in which the states 
that caiTy most of the oscillator strength are clustered about 
~ 850 nm (1.46 eV). Another important difference between 
the two BChl rings is that while the B800s are suiTounded by 
mostly hydrophilic protein residues the binding pocket of the 
B850s is predominantly hydrophobic.^ Thus, although both 
B800s and B850s are chemically identical BChl-a molecules 
their specific spatial arrangement and the nature of their im- 
mediate protein surrounding alter differently their spectral and 
optical properties. For example, it is quite surprising that the 
two peaks, due to the B800 and B850 BChls, in the experi- 
mental OD spectrum of LH2 from Rs. molischianum at room 
temperature ^^-^^ have comparable widths although, as men- 
tioned above, the B800 levels are almost degenerate while the 
B850 levels form a ~ 0.2 eV wide excitonic band. 

Clearly, novel methods for calculating optical spectra of 
PPCs by using computer simulations based entirely on the 
atomic structure of the system would not only provide a better 
understanding and interpretation of the existing experimental 
results but would also help in predicting and designing new 
experiments. The standard procedure to simulate the exper- 
imental spectra of LH2 systems (and PPCs in general) con- 
sists of two step a i - ii i -i - - 7 i. First, the excitation energy spectrum 
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is determined based on the static crystal structure of the sys- 
tem and, second, the corresponding stick spectrum is "dressed 
up" with simulated Gaussian (in case of static disorder) and/or 
Lorentzian (in case of dynamic disorder) line widths char- 
acterized by empirically (and often self consistently) deter- 
mined parameters. Alternatively, the spectral broadening of 
the stick spectrum can also be described through the coupling 
of the electronic excitations (modeled as two- or multilevel- 
systems) to a stochastic heat bath characterized by a model 
spectral density with empirical parameters. While either ap- 
proach may yield excellent agreement between the simulated 
and experimental spectra, the empirical nature of the model 
parameters restricts their predictive power. 

Our proposed method for calculating optical spectra is 
based on a combination of all atom molecular dynamics (MD) 
simulations, quantum chemistry (QC) calculations, and quan- 
tum many-body theory. The conformational dynamics of the 
LH2 ring embedded into its natural environment (a fully sol- 
vated lipid bilayer) are followed by means of classical MD 
simulations. Next, for each BChl, modeled as a quantum 
two level system, the Q,. excitation energy gap and transition 
dipole moment time series are determined along a properly 
chosen segment of the MD trajectory by means of QC cal- 
culations. Finally, the OD and CD spectra are determined 
as weighted sums of the Fourier transform of the quantum 
dipole-dipole correlation function (i.e., the absorption line- 
shape function) which, within the cumulant approximation, 
can be calculated from the sole knowledge of the energy gap 
time series. Formally, this method can also be regarded as a 
two step procedure. First, a stick spectrum is generated from 
the average values of the energy gap time series and, second, 
spectral broadening is applied through the corresponding line- 
shape function weighted by the mean transition dipole (rota- 
tional) strength in the case of OD (CD) spectrum. Since both 
the peak position and the broadening of the optical spectrum 
are obtained from the same energy gap time series determined 
from combined MD/QC calculations, the proposed method re- 
quires no empirical fitting parameters making it ideal for pre- 
dicting optical spectra for PPCs with known structure. Similar 
MD/QC methods were used previously by Mercer et alj^ for 
calculating the OD spectrum of BChl-a in methanol, and by 
Damjanovic et al.— to determine the OD and CD spectra of 
B850s in LH2 from Rs. molischianum. The relationship be- 
tween these studies and the present one will be established 
below. 

The reminder of the paper is organized as follows. The 
theoretical background of the proposed method for calculat- 
ing optical spectra of PPCs is presented in Sec. HI] The em- 
ployed MC simulations and QC calculations are described in 
Sec.|in| The obtained results and their discussion is contained 
in Sec. lIVI Finally, Sec.|y]is reserved for conclusions. 



II. THEORY 

In order to calculate the Unear optical absorption of a PPC 
we assume that the electronic properties of individual pigment 
molecules can be described in terms of a two-level system. 



formed by the ground state and the lowest excited singlet state 
(e.g., the Qy state in the case of BChl-a) involved in the opti- 
cal absorption process. Neglecting for the moment the direct 
interaction between the pigments (e.g., by assuming a suffi- 
ciently large spatial separation between them as in the case of 
the B800s in LH2), we denote these two states for the pig- 
ment (n — 1, . . .A^) as |0) = |0„) and \n) = |1„), respectively. 
Once the interaction between the pigment and its environment 
(composed of protein matrix, lipid membrane and solvent 
molecules) is taken into account these two levels turn into, still 
well separated, energy bands |0;Ao) = |0)|A()) and |n;A,„) = 
|n)|A„), where the quantum numbers Ao and A„ specify the 
state of the n''' pigment on the ground- and excited-state po- 
tential energy surface, respectively. Because the exact quan- 
tum mechanical treatment of the eigenstates |0;Ao)^ \n;Xn) 
and of the corresponding energy eigenvalues S'q^x^^, <Ci,A„ is 
not feasible, usually the quantum numbers Aq and Xn are asso- 
ciated with the vibronic states of the PPC that can be treated 
within the harmonic approximation. Here we follow a differ- 
ent approach in which the dynamics of the nuclear degrees 
of freedom of the PPC are described by all-atom MD simula- 
tions, and the energy gap time series AE„{t) = S'„{t) — (^ci(f) 
is calculated at each MD time step by QC calculations as de- 
scribed below. The main assumption of this approach is that 
the obtained energy gap time series AE„{t) can be used to cal- 
culate approximately equilibrium quantities (such as energy 
gap density of states and time autocorrelation functions) of 
the original system without the knowledge of the exact energy 

gap spectrum Ac^nXM = '^".l,, ^ <^o.k>- 

In the absence of the excitonic coupling between the pig- 
ment molecules, the Hamiltonian of the system can be written 

as J/f = Hq+H, where 

//o = ^|0;Ao)4.;i<,(0;Ao|, (la) 

and 

H = Y,H„=Y,\n-X„KaM,Xn\. (lb) 

The dipole moment operator through which the incident light 
field couples to the n''' pigment is given by 

A« = L <*"A..Ao|n;^)(0;Ao|, (2a) 

where the transition dipole moment (TDM) matrix element 
dn.A„,Afl in Condon approximation'* can be written 

d„.A„,Ao«d„(A„|Ao). (2b) 

Here d„ = (1 |/i„|0} is the real TDM vector whose time series 
can be determined from the same combined MD/QC calcu- 
lations as AE„{t). Note that while (1|0) =0, in general the 
Franck-Condon factors (A„|Ao) are finite"*. 

When the size of the PPC is much smaller than the wave- 
length of the light field, in leading approximation the latter can 
be regarded as homogeneous throughout the system and, ac- 
cording to standard linear response theory, the corresponding 
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OD spectrum is proportional to the dipole-dipole correlation 
function 



(3) 



where fi„.i{t) = e fl„j{Q)e'^°' is the / G {xj^jz} component 
of the time dependent electric dipole operator, and (...) = 
Tr |Zq"' exp(— /3//o) . . .} with )3 = l/kgT the usual temper- 
ature factor and Zq the corresponding partition function. To 
simplify notation, throughout this paper we use units in which 
Ti—\, and apply the convention of implicit summation over 
repeated vector indices. By employing Eqs. Q-Q, after 
some algebra, the quantum dipole correlation function in 
Eq. Q can be expressed as 

Alj(0)A„,/(f )) = d„jd„,,j5„,n {e'"»'e-'""') , (4) 

where 5„„, is the Kronecker delta. By inserting Eq. (0} into 
Eq. one obtains the sought OD spectrum of an aggregate 
of noninteracting pigments in their native environment 



I{co) oc coY,dlA„{(o), 

n 

where the line shape function is defined as 



A„{(0) = Re / dte"^' {e'"<'' e-'""') 



(5a) 



(5b) 



The main difficulty in calculating the quantum time corre- 
lation function in Eq. i5h\ is due to the fact that the Hamilto- 
nians Ho and H„ do not commute. If they would, then the line- 
shape function could be expressed in terms of the energy gap 
density of states (DOS). Indeed, in this case (^e'^^'e^'""') « 
(exp(— /A//„f)), with AH„ = H„ — Hq, and by calculating the 
time integral in Eq. i5h\ would follow 

A„ (co) K 7Z^ (co) , (6a) 
.yK{(o) EE {5{co~AH„)) « {5{co-AE„{t))), (6b) 

where the density of states ,yy{co) is approximated by the 
binned histogram of the energy gap fluctuations AE„{t) ob- 
tained from combined MD/QC calculations. In general, 
Eqs. (|6} overestimate the broadening of the lineshape func- 
tion. Indeed, the Fourier transform of the exact spectral repre- 
sentation of the correlation function 



(7a) 



where p^^ = Zg ^exp{—fiS'Qj^) is the statistical matrix of the 
electronic ground state, yields 

A{co)^2n £ p;i^|(Ao|A„)|25(a)-A4,;L„,Ao) > (7b) 

which can be regarded as a Franck-Condon weighted and ther- 
mally averaged density of state^^. By setting the Franck- 
Condon factors (Ao|A„) equal to unity in i7h\ one obtains 



Eqs. (|6}. Since it is not possible to determine all these fac- 
tors, it is often convenient to use Eqs. (|5Jl as a rough estimate 
of An{(o) for calculating the OD spectrum. 

A systematic way of calculating the correlation function in 
i5h\ is the cumulant expansion method. Here we employ the 
second order cumulant approximation that is often used in op- 
tical spectra calculations^. We have 



(e'^o'e-'^"') = (Texp 



-/ / dTAHjr) 
Jo 



■ exp 



'i{AH„)t- / dT{t-T)%{T) 

Jo 



(8) 



where T is the time ordering operator, AHn{t) = 
e'"'"AH„e-'"o', %,(t) = (dH„{t)dH„{0)}, and dH„{t) = 
AH„{t) — {AH„). To make progress, the quantum statistical 
averages in Eq. (|8j will be approximated with classical ones 
involving the energy gap time series AE„{t), i.e., 

(Af^„) « (A£„(f)) = «„, (9a) 



RQ[^n{t)] « Cn{t) = {5E„{t)5E„{Q)) , (9b) 

where 5E„{t) = AE„{t) — {AE„). While approximating a 
quantum time correlation function by identifying its real 
part with the corresponding classical correlation function 
as in Eq. (|9b} is widely used22i^2i^ other approximation 
schemes have also been considered in the literature^^. Next, 
by invoking the fluctuation dissipation theorem '^,(—0)) = 
exp(-)3a))'^„(a)), where "^(o)) = j'^„dt^„{t)e,x^{im) is 
the Fourier transform of 'rf„{t), the quantum correlation func- 
tion in terms of the real spectral density 



M(0) = \ [in{(0)~%,{-(0)\ = i (X-e-l^") ^n{(0) 



can be written as 



'r„(f) = <(0-<(f) 



(10) 



(11) 



= / — 7„(a))[coth(j3a)/2)cos(Bf-/sin(Bf] . 
Jo It 

By identifying the real part of Eq. (II 1> with Eq. ( I9bt one can 
determine both the spectral density and the imaginary part of 
the quantum correlation function, i.e.. 



7„((b) 2tanh()3a)/2) / dtC„{t) cos m, (12) 
Jo 



and 



K'it)^ [ — Jn((o) sin m. (13) 

Jo 71 

Thus, the lineshape function within the second cumulant ap- 
proximation is 

A„(a)) EE A„(a)-a)„) / e"*'''^ cos[(a) - a)„)f + ^„(f)], 
Jo 

(14a) 
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where the broadening and frequency shift functions are given 
by 

0„(f)= / dT{t-T)C„{T), (14b) 



and 



d(OJn{(0) 



(Ot — Sin Of 

0)2 



(14c) 



A straightforward extension of the above method for calcu- 
lating the lineshape function and the OD spectrum of ex- 
citonically coupled pigment molecules would require the de- 
termination of the energies S'j Xj and TDMs dj correspond- 
ing to the exci tonic states \J;Xj),J= 1,...,N. Unfortunately, 
the required QC calculations (by considering all pigments 
as a single quantum system) are still prohibitively expen- 
sive computationally. Therefore, we employed an effective 
Hamiltonian approximation for determining the time series 
A£'y(f ) = (fy(f) - (?o(f) and djit) from A£'„(f) and d„(f ) of the 
individual pigments. Assuming that these are coupled through 
the usual point dipole-dipole interaction 
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(15) 



In term of the coefficients c^''' = {J\n) the exci tonic TDMs are 



dy = ^(y|n)d„ 



(17) 



Next, by rewriting the Hamiltonian ([^ in diagonal form 
(i.e., in terms of noninteracting excitons) H — Y^jHj — 
\J'^^j)'^j,Xj{J'->^j\' ^fter some algebra one arrives at the 
equations 



(18a) 



L (A,!;.,(0)A„,,(0) = Jldjjdjj {e^"»'e-^"^') . (18b) 



and 



Inserting Eq. ( I18b> into Eq. (|3} one obtains the desired OD 
spectrum of the excitonic system 



where 



Aj{(0) = ReJ^dte"" {e'"<^' e''"'') . 



(19) 



(20) 



We note that by replacing the site index n with the excitonic 
index J most of the above results for noninteracting pigments 
remain formally valid for the corresponding excitonic system 
as well. For example, similar expressions to (|6j and MA\ can 
be easily derived for estimating Aj{(o). 

To conclude this section we derive an expression for the CD 
spectrum of the PPC. By definition, the CD spectrum Icd{(o) 
is the difference between Il{co) and Ir{co), the OD spectra for 
left and right circularly polarized light, respectively. Unlike in 
the case of the OD spectrum, the calculation of Icd{(o) even 
within the leading order approximation requires taking into 
account the spatial variation of the light field across the PPC as 
well as the excitonic coupling between the pigment molecules 
regardless how small this may be. The sensitivity of the CD 
spectrum to geometrical and local details of the PPC makes 
it a quantity difficult to predict by theoretical modeling. The 
CD spectrum is given by^ 



Icd{(0) = ]rM(0) -Ir{(0)] - (oRe f dte" 
4 Jo 

xEy£'-^(r«)^(AL(o)A«,<(0 



where is the relative dielectric permitivity of the medium, 
r„ is the position vector of pigment n, and r„,„ = r„, — r„, the 
eigenvalue equation one needs to solve at every MD timestep 
is 

YjS-^n^nm + Km) " A£'y5„„,]cJ;f^ = . (16) where 



(21) 



where A is the wavelength of the incident light and £, y^. i s the 
unit antisymmetric tensor of rank 3. Inserting Eq. dlSat into 
J21> and making use of Eq. (I20> . we obtain 



/c£>(«) °= (oY^RjAj{(o) 

J 



(22a) 



Rj=yL <-^I") ■ (d« X dm)] HJ) (22b) 



is the so-called rotational strength of the excitonic state /. 
Note, that in the absence of the excitonic coupling all Rj = 
(because for a given J only one coefficient {J\n) is nonzero) 
and the CD spectrum vanishes. The rotational strength plays 
the same role for the CD spectrum as the TDM strength for 
the OD spectrum. Specifically, 7?/ gives the coupling between 
the TDM of the excitonic state J and the orbital magnetic mo- 
ment of the other excitons. The coupling to the local magnetic 
moment is assumed to be small (Cotton effect) and usually is 
discardedii^. 



III. COMPUTATIONAL METHODS 

In order to apply the results derived in Sec. |ll|for calcu- 
lating the OD and CD spectra of the B800 and B850 BChls 
in a single LH2 ring from Rs. moUschianum first we need to 
determine the time series of the Q, energy gap A£'„(£Af) and 
TDM d„(^Af), ^ = 0, 1, . . . ,M, for all individual BChls. We 
accomplish this in two steps. First, we use all atom MD simu- 
lations to follow the dynamics of the nuclear degrees of free- 
dom by recording snapshots of the atomic coordinates at times 
ff = Mf, and then use QC calculations to compute AZi„ and d„ 
for each snapshot. 
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A. Molecular dynamics simulations 

Here we provide a brief description of the simulated LH2 
ring in its native environment as well as the employed MD 
simulation protocol. A more detailed account or the reported 
MD simulations can be found in Ref. 16. A perfect 8-fold LH2 
ring was constructed starting from the crystal structure (pdb 
code ILGH) of Rs. molischianuirfi . After adding the missing 
hydrogens, the protein system was embedded in a fully sol- 
vated POPC lipid bilayer of hexagonal shape. Finally, a total 
of 16 Cl^ counterions were properly added to ensure elec- 
troneutraUty of the entire system of 87055 atoms. In order to 
reduce the finite-size effects, the hexagonal unit cell (with side 
length ~ 60A, lipid bilayer thickness ~ 42A and two water 
layers of combined thickness ^ 35 A) was replicated in space 
by using periodic boundary conditions. The CHARMM27 
force field parameters for proteins^^'^^ and lipids^^ were used. 
Water molecules were modeled as TIP3P^''. The force field 
parameters for BChls and lycopenes were the ones used in 
Ref.[l3- After energy minimization, the system was subjected 
to a 2 ns long equilibration in the NpT ensemble^^ at normal 
temperature {T — 300 K) and pressure ip = l atm), using peri- 
odic boundary conditions and treating the full long-range elec- 
trostatic interactions by the PME method^^. All MD simula- 
tions were preformed with the program NAMD 2.5^^, with a 
performance of 8.5 days/ns on 24 CPUs of an AMD 1800+ 
Beowulf cluster During equilibration an integration time step 
of 2 fs was employed by using the SHAKE constraint on all 
hydrogen atoms'*". After the 2 ns equilibration a 1 ps produc- 
tion run with 1 fs integration step was carried out with atomic 
coordinates saved every other timestep, resulting in A', = 500 
MD snapshots with Af = 2 fs time separation. These config- 
uration snapshots were used as input for the QC calculations 
described below. 



B. Quantum chemistry calculations 

The time series of the Qy transition energies AEn and dipole 
moments d„ of individual BChls can be determined only ap- 
proximately from the configuration snapshots obtained from 
MD simulations. The level of approximation used is deter- 
mined by: (i) the actual definition of the optically active quan- 
tum system, i.e., the part of the system that is responsible for 
light absorption and needs to be treated quantum mechani- 
cally; (ii) the actual choice of the QC method used in the cal- 
culations; and (iii) the particular way in which the effect of the 
(classical) environment on the quantum system is taken into 
account in the QC calculations. Because the optical properties 
of BChls are determined by the cyclic conjugated ;r-electron 
system of the macrocycle the quantum system was restricted 
to a truncated structure of the BChl-a containing 49 atoms in 
the porphyrin plane. The truncation consisted in removing the 
phytyl tail and in replacing the terminal CH3 and CH2CH3 
groups on the macrocycle with H atoms in order to satisfy 
valence requirements. Similar truncation schemes have been 
employed previously^^-*'. In these studies the phytyl tail of 
the BChls was removed but the number of atoms retained in 
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FIG. 1: Normalized DOS, -^(w), for (a) B800, and (b) B850 BChls 
in LH2 of Rs. Molischianum computed as binned histograms of the 
corresponding Q,. excitation energy time series obtained from com- 
bined MD/QC simulations. Whether the charge fluctuations of the 
BChls' environment are included (solid lines) or not (dashed line) 
makes an important difference in ./K(ft)) only for B800. In (b) the 
DOS of the B850 excitons is shown as a thick solid line. 



the optically active macrocycle was different. For example, 
Cory et al.^' in calculating the excited states of the B800 oc- 
tamer and the B850 hexadecamer of LH2 from Rs. molischi- 
anum used 44 macrocycle atoms, while Mercer et al.— in cal- 
culating the absorption spectrum of BChl-a in ethanol used 84 
atoms. According to the crystal structure^, the B800 and B850 
BChls in LH2 from Rs. molischianum differ only in the length 
of their phytyl chain, having a total of 107 and 140 atoms, 
respectively. The removal of the phytyl tail reduces dramati- 
cally both the size of the quantum system and the correspond- 
ing QC computational time. Furthermore, for the truncated 
BChls the non trivial task of automatic identification of the 
Qy excited state in the case of a large number of such compu- 
tations becomes easier and more precise. Although in general 
the different truncation schemes yield excitation energy time 
series with somewhat different (shifted) mean values, the cor- 
responding energy fluctuations, which play the chief role in 
calculating the optical absorption properties of PPC at room 
temperature in their native environment, are less sensitive to 
the actual size of the truncated pigment. 

The Qy excitations of the truncated BChls were calcu- 
lated by using Zerner's semiempirical intermediate neglect 
of differential overlap method parametrized for spectroscopy 
(ZINDO/S) within the single-point configuration interaction 
singles (CIS) approximation^^-*^. Because it is much faster 
and more accurate than most of the computationally afford- 
able ab initio QC methods (e.g., the Hartree-Fock (HF) CIS 
method with the minimal STO-3G* basis set), ZINDO/S 
CIS has been extensively used in the literature to com- 
pute low lying optically allowed excited states of pigment 
moleculesiMMMl. In cases like ours, where thousands of QC 
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FIG. 2: Average transition dipole moments {dj) corresponding to the 
J = 1, . . . , 16 B850 excitonic states. Both (dj) and the correspond- 
ing error bars are expressed relative to the mean dipole moment of 
individual B850s. 



calculations are required, the proper balancing between speed 
and accuracy is absolutely essential. To further increase the 
computational speed, the active orbital space for the CIS cal- 
culations was restricted to the ten highest occupied (HOMO) 
and the ten lowest unoccupied (LUMO) molecular orbitals. 
According to previous studies^^, as well as our own testings, 
the choice of a larger active space has negligible effect on the 
computed states. Our ZINDO/S calculations were carried 
out with the QC program packages HyperChem/& and GAUS- 
SIAN 98^^. In each calculation only the lowest four excited 
states were determined. Only in a small fraction (< 5%) of 
cases was the indentification of the Qj excited state (charac- 
terized by the largest oscillator strength and corresponding 
to transitions HOMO^LUMO and HOMO-1^LUMOh-1) 
problematic requiring careful inspection. We have found that 
even in such cases the Q,. state had the largest projection of the 
TDM along the y-axis, determined by the NB and ND nitrogen 
atoms. 

The effect of the environment on the quantum system was 
taken into account through the electric field created by the par- 
tial point charges of the environment atoms, including those 
BChl atoms that were removed during the truncation process. 
Thus, the dynamics of the nuclear degrees of freedom (de- 
scribed by MD simulation) have a two-fold effect on the fluc- 
tuations of the Qy state, namely they lead to: (1) conforma- 
tional fluctuation of the (truncated) BChls, and (2) a fluctu- 
ating electric field created by the thermal motion of the cor- 
responding atomic partial charges. In order to assess the rel- 
ative importance of these two effects the time series AE„{t) 
were calculated both in the the presence and in the absence 
of the point charges. Since the ZINDO/S implementation in 
GAUSSIAN 98 does not work in the presence of external 
point charges, these calculations were done with HyperChem. 
The ZINDO/S calculations without point charges were car- 
ried out with both QC programs and yielded essentially the 
same result. For each case, we have performed a total of 
12,000 (500 snapshots x 24 BChls) ZINDO/S calculations. 
On a workstation with dual 3GHz Xeon EM64T CPU it took 
~ 2.3 min/CPU for each calculation with point charges, and 
only ~ 0.7 min/CPU without point charges. Thus, on a clus- 
ter of five such workstations, all 24,000 ZINDO/S runs were 
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FIG. 3: Absorption spectrum 7/505(0)) of LH2 foi Rs. molischianum 
calculated as a combined DOS of B800 BChls and B850 excitons 
weighted by the corresponding dipole strengths (solid line). /Dos(ft)) 
was blueshifted by 20 meV in order to overlay its B850 peak with the 
corresponding one in the experimental OD spectrum^^ (dashed line). 



completed in ^ 1 .9 + 0.6 = 2.5 days. 



IV. RESULTS AND DISCUSSION 

The time series of the Qy excitation energies A£„(f(') and 
TDMs d„(f£), (f£ ^ Mf; 1^0,... ,N,; N, ^ 499; Af = 2 fs), 
were computed with the ZINDO/S CIS method, described 
in Sec. HiTbI for both B850 (n = 1 , . . . , 16 ) and B800 (« = 
17, . . . ,24) BChls in a LH2 ring from Rs. molischianum, us- 
ing snapshots from the all atom MD simulation described in 
Sec, nil AI The calculations were done both with and with- 
out the point charges of the atoms surrounding the truncated 
BChls. The 1 ps long time series appear to be sufficiently 
long for calculating the DOS of the Qy excitation energies and 
the coiTesponding OD and CD spectra. In Ref. 28 it has been 
found that at least a 2.2 ps long MD trajectory was needed 
for proper evaluation of optical observables related to their 
MD/QC calculations. However, our test calculations showed 
no significant difference between the energy gap autocorrela- 
tion functions calculated from a 2 ps and a 1 ps long energy 
gap time series. Therefore, to reduce the computational time 
we have opted for the latter 



A. Energy gap density of states (DOS) and transition dipole 
moments (TDMs) 

FigureQlshows the Qy energy gap DOS, ,yV{(i)), of the in- 
dividual B800 [top (a) panel] and B850 [bottom (b) panel] 
BChls calculated, according to Eqs. (|6}, as normalized binned 
histograms of the time series AEgsoo = ^nik) with n = 
17, . . . ,24, and AEbssq = AE„{t() with n = 1, . . . , 16, respec- 
tively. In order to eliminate the noise due to finite sam- 
pling, the graphs have been smoothened out by a running 
average procedure. The same smoothing out procedure has 
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FIG. 4: Normalized autocorrelation function C{t) /C(0) of the energy 
gap fluctuations 8E{t) = E{t) — (E) for individual B800 (dashed 
line) and B850 (solid line) BChls, calculated using Eq. . The mean 
square energy gap fluctuations are 

Cbsoq{0) = 3.16 X 10-3 eV2 and 

Ci}85o(0) = 8.68 X lO-"* eV2 



been applied to all subsequent lineshape and spectra calcula- 
tions. In the absence of the point charge distribution of the 
environment ,yy{co) for B800 and B850 (dashed lines) are al- 
most identical, having peak position at 1.51 eV (817 nm) and 
1.515 eV (818 nm), and full width at half maximum (FWHM) 
5 1 meV and 59 meV, respectively. The fact that the peak po- 
sition practically coincides with the mean energy gap is in- 
dicative that the DOS is symmetric with respect to its maxi- 
mum. It should be noted that essentially the same mean en- 
ergy gap of 1.5 eV was obtained in similar MD/QC calcu- 
lations (i) by us (data not shown) in the case of a truncated 
BChl-a in vacuum but artificially coupled to a Langevin heat 
bath at room temperature, and (ii) by Mercel et al.^^ for a 
BChl-a solvated in methanol also at room temperature. Al- 
though in case (ii) the width of the DOS appears to be some- 
what broader (FWHM« 65 eV) than in case (i), for which 
FWHM« 58 meV, based on these results one can safely con- 
clude that the thermal motion of the nuclei in individual BChls 
lead to Qy energy gap fluctuations that are insensitive to the 
actual nature of the nonpolar environment. Since in LH2 
from Rs. molischianum the surrounding of the B800s is po- 
lar while that of the B850s is not, one expects that once the 
point charges of the environment are taken into account in the 
QC calculations J/{(£)) should change dramatically only in 
the case of B800. Indeed, as shown in Fig. ^ (solid line), 
in the presence of the point charges the peak of o/(^85o(g)) is 
only slightly red shifted to 1.502 eV (825 nm) and essentially 
without any change in shape with FWHM« 53 meV. By con- 
trast, the DOS for B800 in the presence of the point charges 
(Fig. nj) has quahtatively changed. The induced higher en- 
ergy fluctuations not only spoil the symmetry of o/>^8oo(ro) 
but also increase dramatically its broadening, characterized 
by FWHM« 100 meV. Thus, in spite of a small blueshift to 
1.528 eV (811 nm) of the peak of ./^B8oo(fi^) the mean value 
of the energy gap (AZifigoo) = 1.556 eV (797 nm) is increased 
considerably, matching rather well the experimental value of 
800 nm. 
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FIG. 5: Spectral density function ^(h') for B800 (dashed line) and 
B850 (solid line) obtained according to Eq. <12t . 



The time series of the excitonic energies A£y(rf), J ~ 
1 , . . . , 16, of the B850 BChls were determined by solving for 
each MD snapshot, within the point-dipole approximation, the 
eigenvalue equation (I16> . In calculating the matrix elements 
( I15> r„ was identified with the position vector of the Mg atom 
in the «-th BChl. Consistent with the Condon approximation, 
the magnitude of the computed B850 TDM time series exhib- 
ited a standard deviation of less than 4% about the average 
value (t/flgso) = 11 .77 D. The latter is by a factor ofk= 1 .87 
larger than the experimentally accepted 6.3 D value of the Qy 
TDM of BChl-a^^. Thus, to account for this overestimate of 
the TDM by the ZINDO/S CIS method, instead of using a 
reasonable value of 1.86 for the relative dielectric constant of 
the protein environment in Eq. (I15> we set for all the calcula- 
tions reported in this paper = 6.5 (= 1.86 x 1.87^). Equiv- 
alently, one can rescale all TDMs from the ZINDO/S calcu- 
lations by the factor fc"' and set — 1.86. Either way the 
mean value of the nearest neighbor dipolar coupling energies 
between B850s were 27 meV« 220 cm^' within a protomer 
and 24 meV« 196 cm^' between adjacent heterodimers. Just 
like in the case of individual BChls, the DOS coiTesponding 
to the B850 excitonic energies (Fig.^ - thick line) was cal- 
culated as a binned histogram of AEj{t(). As expected, the 
excitonic DOS is not sensitive to whether the point charges of 
the environment are included or not in the B850 site energy 
calculations. 

The mean excitonic TDMs, calculated from Eq. and 
expressed in terms of (dBSSo), are shown in Fig.|2] The er- 
ror bars represent the standard deviation of the time series 
dj{t(). In agreement with previous studies, most of the dipole 
strength is amassed into the lowest three excitonic states. 

As discussed in Sec. HH a rough estimate of the lineshape 
function can be obtained as the combined DOS of the B800 
BChls and B850 excitons. In this approximation the OD spec- 
trum reads 
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FIG. 6: Lineshape functions AB8()o(Aa)) (dashed line) and 
Ag85()(Aa)) (solid line). 



FIG. 7: Computed (solid line) and experimental (dashed line) absorp- 
tion spectra (in arbitrary units) of the BChl aggregate in Rs. MoUschi- 
anum LH2. The computed spectrum has been blue shifted by 20 meV 
for best match. 



where the B800 index in the last term means summation over 
all B800 BChls. Figure |3] shows the calculated Idos{(») 
blueshifted by 20 eV (solid line) in order to match the 
B850 peak position with the one in the experimental OD 
spectrum. '^-^^ (dashed line). While the B850 band and the rel- 
ative heights of the two peaks in Idos{(o) match rather well 
the experimental data, the position and the broadening of the 
B800 peak do not. This result clearly shows that in general 
peak positions in optical spectra may be shifted from the cor- 
responding peak positions in the excitation energy spectrum 
due to correlation effects between the ground and optically 
active excited states. The latter may also lead to different line 
broadening of the corresponding peaks. Thus, it appears that 
in principle, methods for simulating optical spectra in which 
the position of the peaks are identified with the computed ex- 
citation energies (stick spectrum) are not entirely correct and 
using instead more sophisticated methods that include quan- 
tum correlation effects should be preferred. Such method, 
based on the cumulant approximation of the lineshape func- 
tion as described in Sec. HI] is used in the next section for 
calculating the OD spectrum of an LH2 ring from Rs. molis- 
chianum. 



B800 [B850] BChls according to the formula 



B. Absorption (OD) spectrum 



The key quantity for calculating the lineshape functions of 
the individual B850 and B800 BChls is the (classical) auto- 
coiTelation function C„(f) = {dE„{t)dE„{0)) of the energy 
gap fluctuation 5£'„(f) — AE„{t) — {AE„) determined from 
the combined MD/QC calculations. Because the time series 
AE„{t) were too short for a proper evaluation of the ensemble 
average in the individual C„(f ), a single time correlation func- 
tion CB8oo(f ) [CB&5o{t)] was determined by averaging over all 



1 



where M 
andM 



1 N,-e 
- — - 5E„{te + tk)5E,„{tk) 



8, m== 17,. 
: 16, m = 1, . 



,24 for a =B800, 
,16 for a =B850. 



(24) 



The normalized correlation functions Ca{t) /Ca{Q), OL G 
{B800,B850}, are plotted in Fig.H Ca(0) = (Sfi^) represents 
the variance of the energy gap fluctuations with Cb80o(0) = 
3.16 X 10-^ eV2 and Cb85o(0) = 8.68 x 10-"* eV^. The be- 
havior of the two correlation functions is rather similar during 
the first 150 fs. Following a sharp decay to negative values 
in the first 9 fs, both functions exhibit an oscillatory compo- 
nent of approximately 18.5 fs period and uneven amplitudes 
that, in general, are larger for the B800. After ^ 150 fs, the 
autocorrelation functions behave in a distinctive manner, both 
becoming negligibly small for t > 400 fs. 

The spectral densities Ja{(o) for B800 and B850, deter- 
mined according to Eq. (I12> . are shown in Fig.|5l The promi- 
nent peak about (o,, ~ 0.22 eV is due to the fast initial decay 
of Ca{t)- Being reported in previous studies'^'^^, by using 
both ab initio (HF/CIS with STO-3G* basis set) and semi em- 
pirical QC methods, these spectral features appear to be in- 
trinsic properties of BChl-a, most likely originating from a 
strong coupling of the pigment to an intramolecular C=0 vi- 
bronic mode. Often, the environment in a PPC is modeled as 
an equivalent harmonic (phonon) heat bath for which the cu- 
mulant approximation is exact^. The corresponding phonon 
spectral density can be written asJ{(o) — ay^'^xgj^5{(0 — (Oi), 
where gi is the coupling constant to the phonon mode X. 
Thus, one can interpret the magnitude of the spectral functions 
in Fig.|5]as a measure of the coupling strength to phonons of 
that particular frequency. The complex structure of the spec- 
tral functions indicate that all inter and intra molecular vi- 
bronic modes with frequency below cOp will contribute to the 
lineshape function. Hence, attempts to use simplified model 
spectral functions appear to be unrealistic even if these may 
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FIG. 8: Mean rotational strength of the excitonically coupled B800 
(circles) and B850 (rectangles) BChls as a function of the corre- 
sponding excitonic energies. The purpose of the thin lines are to 
guide the eye. 



lead to absorption spectra that match the experimental results. 

The lineshape functions of individual B800 and B850, cal- 
culated from Eqs. il4\ . are plotted in Fig.|6] The origin of the 
frequency axis corresponds to the mean energy gaps CObhoo and 
0^850. respectively. The highly polarized surrounding of the 
B800 BChls in Rs. molischianum renders AB8oo(ft)) twice as 
broad (FWHM« 26 meV) as Ab%sq{(o) (FWHM« 13 meV). 
Also, the redshift of the peak of the former (Ao w 25 meV) 
is more than three times larger than that of the latter (Ao) w 
7 meV). 

Since the available simulation data is not sufficient to prop- 
erly estimate the excitonic lineshape functions Ay (©), by ne- 
glecting the effect of exchange narrowing^'^-^'' , we approxi- 
mated these with AB%5i){(£>). Thus, the OD spectrum of the 
LH2 BChls was calculated by using the formula 



/(co) o= O) 



COj) + 8c/|800^B800((» - C0b800) 



(25) 

where (Oj = (AEj). 

As shown in Fig. after an overall blueshift of 20 meV, 
I{co) matches remarkably well the experimental OD spectrum, 
especially if we take into account that it was obtained from the 
sole knowledge of the high resolution crystal structure of LH2 
from Rs. molischianum. The reason why both B800 and B850 
peaks of /(o)) are somewhat narrower than the experimental 
ones is most likely due to the fact that the effect of static dis- 
order is ignored in the present study. Indeed, our calculations 
were based on a single LH2 ring, while the experimental data 
is averaged over a large number of such rings. While compu- 
tationally expensive, in principle, the effect of static disorder 
could be taken into account by repeating the above calcula- 
tions for different initial configurations of the LH2 ring and 
then averaging the corresponding OD spectra. 

To conclude this section we would like to relate the present 
work to previous two combined MD/QC studies In Mer- 
cer et al.^^ it is argued that the ab initio QC method (HF/CIS 
with the STO-3G* basis set) should be preferred to semi em- 
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FIG. 9: (a) CD spectrum contributions due to B800 (dashed line) 
and B850 (solid line) BChls. (b) Comparison between the computed 
(solid line) and experimental CD spectrum of the BChl aggregate in 
Rs. Molischianum LH2. 



pirical methods for calculating optical spectra because it re- 
produces better their experimental results. The FWHM of 
their calculated semi empirical and ab initio absorption spec- 
tra of BChl-a in methanol are ^ 65 meV and 125 meV, 
respectively. These values are similar to the ones we obtained 
for the same type of calculations for BChl-a in vacuum and 
in LH2 embedded in its native environment. Since except 
Ref. I23 all experimental results on the absorption band 
of BChls we are aware of have a FWHM of < 80 meV at 
room temperature, we conclude that in fact the semi empiri- 
cal ZINDO/S method should be preferable to the ab initio QC 
method. In general, the latter overestimate the broadening of 
the OD spectrum by a factor of 2 to 3. In Damjanovic et al. ^4, 
the OD spectrum of individual B850s (i.e., without excitonic 
coupling) calculated with the ab initio method yielded the 
same FWHM of ^ 125 meV as in Mercer et al.^**. Once 
the excitonic coupling was included within the framework of 
a polaron model, and it was assumed that the entire oscilla- 
tor strength was carried by a single exciton level of a perfect 
B850 ring, the FWHM of the resulting OD spectrum was re- 
duced to ^ 43 meV as a result of exchange narrowing. Thus 
the obtained OD spectrum for the B850s appeared to match 
rather well the corresponding part of the experimental one. 
However, by applying the same method to the B800s, where 
there is no exchange narrowing, the polar environment fur- 
ther broadens the corresponding OD spectrum to a FWHM of 
^ 250 meV that is clearly unphysically large. Thus, the con- 
clusion again is that the ZINDO/S CIS semi empirical method 
should be preferred for calculating optical spectra of PPCs. 
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C. Circular dichroism (CD) spectrum 

The CD spectrum of the LH2 BChls from Rs. molischi- 
anum was determined by following the theoretical approach 
described at the end of Sec. HI] and by employing the same 
time series (obtained from the combined MD/QC calcula- 
tions) used for calculating the OD spectrum. 

First, the rotational strength of both B850 excitons and 
B800 BChls were determined by using Eq. M2h\ . In this equa- 
tion, just like in the case of the point-dipole interaction matrix 
elements M5\ . the vector r„ described the position of the Mg 
atom in the n''^ BChl. As already clarified in Sec.|lllthe cal- 
culation of the rotational strength of the B800 BChls requires 
solving the corresponding excitonic Hamiltonian M6\ regard- 
less of how small the dipole-dipole coupling is between these 
BChls. The calculation does not yield either noticeable cor- 
rections to the B800 excitation energies or admixture of the 
corresponding Q, states, however, it leads to sizable mean ro- 
tational strengths as shown in Fig. |S] (filled circles). Similarly 
to the TDM strengths (Fig.|2}, the largest (negative) mean ro- 
tational strengths are carried by the four lowest B850 exci- 
tonic states as shown in Fig. |8l (open squares). The second 
highest excitonic state also has a sizable rotational strength 
and is responsible for enhancing the positive peak of the B800 
contribution to the CD spectrum (Fig.|9^). 

Second, the CD spectrum is calculated from Eq. (|Sa) 
where the summation index J runs over all B850 and 
B800 excitonic states and Ay(co) = Aa{(0 — coj), with a e 
{B850,B800}. Figure |9ji shows the CD spectrum contribu- 
tion by the B850 (solid line) and B800 (dashed line). Both 
contributions have the same qualitative structure with increas- 
ing energy: a pronounced negative peak followed by a smaller 
positive one. The B850 negative CD peak is about twice as 
large as the corresponding B800 peak. The total CD spec- 
trum, given by the superposition of the B850 and B800 con- 
tributions, is shown in Fig.|9j) (solid line) and matches fairly 
well the experimental spectrum.'^ (dashed line). It should be 
emphasized that apart from an overall scaling factor the CD 
spectrum was calculated from the same MD/QC data as the 
OD spectrum by following the procedure described above. 



the LH2 ring from Rs. molischianum the OD and CD spec- 
tra of this PPC can be predicted with reasonable accuracy 
at affordable computational costs. The configuration snap- 
shots taken with femtosecond frequency during the MD sim- 
ulation of the PPC in its native, fully solvated lipid-membrane 
environment at room temperature and normal pressure pro- 
vide the necessary input for the QC calculations of the opti- 
cal excitation energies and transition dipole moments of the 
pigment molecules. The obtained time series are used to 
evaluate within the second cumulant approximation the op- 
tical lineshape functions as the Fourier transform of the quan- 
tum dipole-dipole correlation function. Our choice of the 
ZINDO/S CIS method for the QC calculations was motivated 
by the fact that it is almost two orders of magnitude faster and 
much more accurate than the most affordable ab initio method 
(HF/CIS with the STO-3G* basis set). Compared to the for- 
mer, the latter method overestimates by a factor of 2 to 3 both 
the excitation energies and the broadening of the energy spec- 
trum. Just like in several previous studies we have 
found that the ZlNDO/S method repeatedly yields results in 
good agreement with existing experimental data. 

By investigating the excitation energy spectrum of the LH2 
BChls both in the presence and in the absence of the atomic 
partial charges of their environment we have convincingly 
demonstrated that the large broadening of the B800 peak is 
due primarily to the electric field fluctuations created by the 
polar surrounding environment of the B800s. There is no 
such effect for the B850s which sit in a nonpolar local en- 
vironment. The broadening of the B850 peak is due to the 
sizable excitonic coupling between these BChls. Since only 
the lowest three excitonic states carry most of the available 
dipole strength, in spite of the ^ 0.2 eV wide excitonic band, 
the B850 absorption peak has a FWHM only slightly larger 
than the B800 one. 

It is rather remarkable that both the OD and the CD spec- 
tra of the considered LH2 complex are fairly well predicted 
by our combined MD/QC method. However, a more thor- 
ough testing of the proposed method, involving other PPCs is 
necessary to fully establish its capability of predicting optical 
properties by using only atomic structure information. 



V. CONCLUSIONS 

The continuous increase in processor power and availabil- 
ity of high performance computer clusters, along with the 
growing number of high resolution crystal structures of mem- 
brane bound pigment-protein complexes make feasible the 
theoretical characterization of the spectral and optical proper- 
ties of such systems at atomic level. By applying an approach 
that combines all atom MD simulations, efficient semi em- 
pirical QC calculations and quantum many-body theory we 
have shown that starting solely from the atomic structure of 
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